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Abstract 

Originally a model for wave propagation on the line, the Toda lattice is 
a wonderful case study in mechanics and symplectic geometry. In Flaschka’s 
variables, it becomes an evolution given by a Lax pair on the vector space 
of real, symmetric, tridiagonal matrices. Its very special asymptotic behavior 
was studied by Moser by introducing norming constants, which play the role 
of discrete inverse variables in analogy to the solution by inverse scattering of 
KdV. It is a completely integrable system on the coadjoint orbit of the upper 
triangular group. Recently, bidiagonal coordinates, which parameterize also 
non-Jacobi tridiagonal matrices, were used to reduce asymptotic questions 
to local theory. Larger phase spaces for the Toda lattice lead to the study of 
isospectral manifolds and different coadjoint orbits. Additionally, the time one 
map of the associated flow is computed by a familiar algorithm in numerical 
linear algebra. 

The text is mostly expositive and self contained, presenting alternative 
formulations of familiar results and applications to numerical analysis. 

Keywords: Toda lattice, completely integrable systems, QR algorithm, isospectral 
manifolds, inverse scattering. 
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1 Introduction 

The Toda lattice is a beautiful mathematical object: mathematical miracles and 
serendipity are abundant. Somehow, everybody has something to say about it, 
and this extends way beyond mathematicians: the connections with physics and 
numerical analysis are fruitful and clarifying. Moreover, the formalism is versatile, 
accomodating a number of interesting dynamical systems. 

In this text, we present some of these many aspects of the Toda lattice. After 
more than thirty years of intriguing discoveries, the subject is not exhausted even 
in its most elementary formulation. We aim at concrete examples, and provide 
occasional pointers to more abstract approaches. 

In Section 2, the differential equation is presented as the physical system intro¬ 
duced by Toda and converted into an evolution on Jacobi matrices by Flaschka’s 
remarkable change of variables. The two basic properties are discussed: the con¬ 
servation of eigenvalues and the very simple asymptotic behavior. A comparison 
with the archetypical Lax pair, the Korteweg-de Vries equation, motivates the in¬ 
troduction of Moser’s inverse variables on Jacobi matrices, given by eigenvalues and 
(discrete) norming constants. 
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author gratefully acknowledges support from CAPES, CNPq and FAPERJ, Brazil. 
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The change of variables taking the original mechanical system in to the set of 
Jacobi matrices with a fixed trace j7o leads to a brief description of the Hamiltonian 
formalism. After a brief review of completely integrable systems and coadjoint orbits 
as symplectic spaces, we interpret the Toda equation as a Hamiltonian system on 
coadjoint orbits (in particular j7o), as done originally by Adler and Kostant. We 
then consider larger phase spaces: generic coadjoint orbits and isospectral manifolds 
of tridiagonal matrices. 

Section 4 is another geometric approach to the Toda lattice, which might have 
anticipated its discovery by decades. It leads very naturally to the solutions by 
factorization by Symes. Section 5 is dedicated to an interpolation theorem which 
connects the subject to numerical analysis: one Hamiltonian in the completely in¬ 
tegrable collection associated to the Toda equation gives rise to a flow which, at 
integer times, consists of matrices obtained by the so called QR algorithm, exten¬ 
sively used in algorithms to compute eigenvalues. 

Norming constants have a drawback: they do not extend to the boundary of the 
set of Jacobi matrices, which is where the fine asymptotic properties of the flows 
occur. Section 6 describes bidiagonal coordinates, which provide charts for the full 
isospectral manifold of tridiagonal matrices. A brief interlude on QR algorithms 
with shifts is provided as an example of the versatility of this new instrument. 

As shown above, the content of this review is biased. Indeed, a short text 
cannot provide coverage for all the ramifications of the subject. As a minimal list 
of alternatives, the curious reader is invited to consider [32 and [35] for algebraic 
aspects related to integrability, m for uses of symplectic geometry to the analysis of 
(variations of) the original Toda system, [23] for an essentially orthogonal overview 
and |37] for a relativistic mutation, one of many interesting physical directions from 
our starting point. 


2 Physical origins and Flaschka’s variables 

In the mid sixties, the Japanese physicist M. Toda [42] proposed a model for wave 
propagation along n particles in a line by the Hamiltonian 

^ n n —1 

Hix,y) = + 

1 1 

As usual, positions Xk and velocities yk vary as 

Xk = = Vk, Vk = = exp(a;fc_i - Xk) - exp(a;fc - Xk+i), i = 

ayk oxk 

Here, xq = —oo, Xn+i = +oo — this is the non-periodic case of the Toda lattice. The 
particles are labeled and from a mathematical viewpoint, there is nothing wrong 
with occasional collisions: they simply pass each other. 

Clearly, the energy H{x(t),y{t)) is a conserved quantity. The invariance under 
translation in configuration space (i.e., the fact that H{x, y) = H(x -\- c, y), for any 
fixed c S R") implies that the center of mass of the system moves uniformly (i.e., 
the linear momentum of the system is preserved). 

Following Flaschka[16]. shift to center of mass coordinates (i.e., work with dif¬ 
ferences Xk — Xk+i), get rid of exponentials and adjust constants, 

ak = -^, k = l,...,n, i k = l,...,n-l 

and the evolutions for x{t) and y{t) become 

= 2(&fc - ), k = l,...,n, = bkiak+i - ak), fc = 1,... ,n - 1 
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which, surprisingly, can be cast as the matrix differential equation 


j = [j, = j[WskJ) — 


{J) 


Here, J is a Jacobi matrix with diagonal entries and principal off-diagonal entries 
bk- A Jacobi matrix J is a real n x n symmetric matrix which is tridiagonal (i.e., 
Jij = 0 if |j — j| > 1) and such that the principal off-diagonal entries (those for 
which \i — j\ = 1) is a positive number. The matrix IlskM is skew-symmetric, with 
the same lower triangular part as M. 

There is nothing wrong in considering equation (J) for arbitrary real symmetric 
matrices (and even non-symmetric matrices). Explicitely, consider the differential 
equation on real symmetric matrices S(t) given by 


s'it) = [s{t),n,kSit)], 3(0) = So. 


(S) 


2.1 Lax pairs, asymptotic behavior 

We present two fundamental facts. The first one is a question about differential 
equations. What kind of differential equation on matrices gives rise to an evolution 
M{t) which preserves the eigenvalues? More precisely, what kind of vector field 
gives rise to a flow of the form M{t) = {P{t))~^M{0)P{t)l Simply take derivatives. 


M'{t) = -(P(<))"^P'(t)(P(t))"^M(0)P(t) -b (P(t))-^M(0)P'(t) 
-(P(t))-ip'(t)M(t) + M(t)(P(t))-ip'(t) = [M(t), (P(t))-ip'(t)]. 


In words, equations given by Lax pairs M'{t) = [M{t),X{t)\ are solved by 
conjugating the initial condition. In particular, if the initial condition is a real, 
symmetric matrix 5'(0), and we require P to be an orthogonal matrix (so that 
symmetry, together with spectrum, is preserved along the orbit S{t)), one should 
consider a Lax pair of the form S'{t) = [S'(t), A(t)], where A{t) is a real, skew 
symmetric matrix (this follows from the standard computation — take the derivative 
at 0 of a curve of orthogonal matrices Q{t) with (5(0) = I). 

Proposition 1 The solution S{t) of (S) starting from a real, symmetric matrix 
S'(O) is of the form Sit) = Q{t)*S{i))Q{t): it is well defined for all t G R. 

An explicit form of Q{t) is given in Section 

Proof: We are left with showing global existence: the matrix norm | [S'! P = tr S'^S 
is conserved along an orbit. 

We now consider the asymptotic behavior for S{t) and J{t), the orbits starting 
from arbitrary symmetric and Jacobi matrices. In the Jacobi case, we may get 
information from the physical interpretation: informally, minimizing the potential 
energy (the second term in the Hamiltonian) leads to the spreading of the particles, 
so that one might expect Xk —Xk+i —oo. Thus, for t —>• cjo, particles should be es¬ 
sentially independent from each other and undertake uniform motion. In Flaschka’s 
variables, the orbit Jit) should converge to a diagonal matrix. 

Proposition 2 The limits fort —^ ±oo of Sit) are diagonal matrices. Whent —> oo 
(resp. t —>■ —oo), J{t) converges to a diagonal matrix with strietly decreasing (resp. 
increasing) entries along the diagonal. 

Proof: The differential equations for the diagonal entries of S{t) are given by 
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S'22 = 2(-^|l + ^2^3 + • ■ • + SL), 

and, in general, 5'^^, = — ^Ij- The equation for indicates that 

it is nondecreasing in time. The same is true not for S '22 but for S'Ji + S 22 , and 
in general the partial traces J2j=i 8 -'^® nondecreasing. Since the norm of S(t) is 
constant, all entries Sij and the partial traces (and their derivatives) are uniformly 
bounded. In particular, the derivatives of the partial traces are integrable, Lipschitz 
functions on R and hence all entries Sij,i 7 ^ j must go to zero. This implies the 
diagonal convergence of S{t) (and hence of J{t)) at ± 00 . 

From equation b'^. — bk{ak+i — Ofe ), the entries 6 fc(t) of J{t) are always positive 
(indeed, if = 0 at some time, it is always zero). When t —^ 00 , we must then have 
Ofc+i — Ofe < 0, so that the diagonal limit matrix has entries in strictly decreasing 
order (the fact that the eigenvalues of a Jacobi matrix are all distinct is proved in 
Lemma IQ) . A similar argument obtains the result for t —>■ — 00 . 

As a side remark, notice that the above proposition implies the spectral theorem 
(only in finite dimensions: the presence of continuous spectrum brings up new issues 
to the asymptotics of the Toda flow my The interested reader might enjoy a proof 
of the Wielandt-Hoffman theorem along similar lines l|13jl. 

2.2 Scattering and inverse variables 

Flaschka was certainly inspired by P. Lax’s casting of the Korteweg-de Vries (KdV) 
equation into Lax pair form. Briefly, the evolution 

Ut(t, X) = 6 u(t, x) Ux{t, x) - Uxxx{t, x), f > 0, cc G K 

is equivalent to the operator evolution 

L'it) = [L{t),Hit)] , 

where 

L{t)f{x) = -f'\x)+u{t,x)f{x), 

H{t)f{x) = ‘if"'{x) - 6u{t,x)f'{x) - 3ux;it,x)f{x). 

Lax then reproved the seminal discovery of Gardner, Greene, Kruskal and Miura 
([S]) that the the evolution of KdV varies the operator L preserving its eigenvalues. 
Analogous computations replacing L by higher order operators lead to the Gelfand- 
Dickey flows, which include the Boussinesq equation and more ( [5] ). 

The connection between differential equations like KdV and the inverse scat¬ 
tering method discovered in the late sixties l|19j. m) led to intense research on 
integrable systems. Originally as a formalism, inverse data provided (infinite di¬ 
mensional) action-angle variables for the KdV equation. 

Inverse variables for Jacobi matrices were pervasive in the early approaches to 
the Toda flows, starting with Moser ([32]). The theorem below, that he attributed 
to Stieltjes, provides inverse variables for the Toda equation (J) ([36], [21]!. 

Lemma 2.1 Jacobi matrices always have distinct real eigenvalues. The first and 
last coordinate of each eigenvector is a nonzero number. 

Proof: Suppose Jv = Xv and z;[l] = 0: equating the first coordinates in the equa¬ 
tion obtains v\f2\ = 0, and successively — thus, the extreme entries of eigenvectors 
are nonzero. Now, two independent eigenvectors associated to the same eigenvalue 
would yield an eigenvector with first coordinate equal to zero. 
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Theorem 2.2 The set of Jacobi matrices is diffeomorphic to the product 

{(Ai, A 2 ,..., A„); Ai < A 2 < ... < A„} X {{vi,V2,... ,Vn) ]v^ > 0^ '^vf = 1}. 

i 

More precisely, the diffeomorphism takes J to its ordered eigenvalues and to the 
first coordinates of its associated eigenvectors, normalized so as to be positive. 

Parlett made the intriguing observation that this algorithm is used by numerical 
analysts to obtain a tridiagonal matrix from a full symmetric matrix, not from a 
diagonal matrix. 

Proof: We sketch a procedure to invert this map presented in m- Define the 
matrix A and the vector v as 

A = diag(Ai, A 2 ,..., A„) , v = {vi,V 2 , ■ ■ ■ ,Vn) 

and consider the sequence of vectors v, Av, ..., A'^~^v. The vectors qi = v,q 2 , ■ ■ ■ ,qn 
obtained from this sequence by applying the Gram-Schmidt method are the columns 
of an orthogonal matrix Q for which J = AQ. 

The entries of the vector v are the norming constants associated to J. Still in 
analogy with KdV, as J{t) solves equation (J), the eigenvalues stay put, as Flaschka 
knew, and Moser showed that the norming constants v(t) varied in a simple fashion: 
simply normalize (under the Euclidean norm) the vector exp(tA) u(0). 

From Proposition [5J at t = ±cxd Jacobi orbits J{t) converge to diagonal matri¬ 
ces with ordered eigenvalues along the diagonal entries. Using norming constants, 
Moser computed the scattering map of the Toda flow. Particles group in pairs with 
the same asymptotic velocity at extreme times: the quantity of interest is the shift, 
the distance between the two straight lines tangent to the asymptotic motion at 
±cxD of particles in the same pair. 

Norming constants, alas, do not parameterize the limit matrices of the flows. 
Actually, they degenerate on a large part of the boundary of the set of Jacobi 
matrices. This will be circumvented by the bidiagonal coordinates in section [6l 

There is a basic issue which has not been handled carefully so far. Equation 
(J), given in terms of skew symmetric matrices A{t), was shown to imply that 
the solution is an orthogonal conjugation of the initial condition J(0) — thus, 
in particular, J[t) is always real, symmetric, and eigenvalues are preserved. But 
why should the evolution preserve the tridiagonal form? This has to happen, if 
Elaschka’s change of variable preserves the physical meaning of the variables. Also, 
this fits with the solution by inverse variables of (J). A more conceptual argument 
showing that J(t) is always a Jacobi matrix will be presented in the next section. 
Yet another argument will come up in Section [4.21 


3 Some symplectic geometry 

Erom its physical description in terms of positions Xi and velocities yi, it is clear 
the Toda lattice admits a Hamiltonian formulation in By Flaschka’s change 

of variables, after removal of the (trivial) evolution of the center of mass, the phase 
space for the differential equation becomes the set of Jacobi matrices with trace 
equal to zero, a cone of dimension 2n — 2, from Theorem l2.21 Somehow, one should 
be able to transfer the original Hamiltonian formulation to the new variables, and 
still proceed with the study of the Toda lattice as a problem in mechanics within 
the new phase. This indicates a more general context, the starting point of a vast 
field, symplectic geometry. In the next subsection, we outline some basic requisites 
for this project. 
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3.1 Complete integrability 

We start with a brief description of a more general definition of Hamiltonian for¬ 
mulation of a vector field. Take a differential equation with phase space M, or 
more precisely, consider the associated vector field Z defined on the tangent bundle 
TM. First equip M with a closed, nondegenerate, 2-form w: the pair {M,uj) is 
a symplectic manifold, necessarily of even dimension, say 2n (excellent sources for 
symplectic geometry are [3], PO]!. 

Each Hamiltonian H : M ^ M. gives rise to a vector field Xh as follows: the 
contraction of oj with Xh should obtain the 1-form dH. Said differently, for every 
vector field Y in M, at each m € M 

u)m{XH{rn),Y{m)) = dH{Y{m)) = Y{H){m), 

where Y{H){m) is the derivative of H along Y at m. The vector field Z admits a 
Hamiltonian formulation \i Z = Xh for an appropriate choice of w and H. 

The simplest example is the standard 2-form in = {{x,y); x,y € R"} given 
by w = A yk- Let us find the vector field Xh = Thk-^^^^k + Sfc ^kdyk 

associated to the Hamiltonian H : R^" ^ R (here, dz is the partial derivative in the 
variable z, and Hz = dzH). For an arbitrary vector field Y = J^k ^kdx^ +Sfe ^kdx^ i 

^{Xh , Y') = 'y ( Qk^k T 'y ^ l^k^k — 'y ^ ^^kUtx^ T Y.bkHy,=dH{Y), 

k k k k 

and, as is well known, Xh = ^Vk 

When is a Hamiltonian vector field Xh completely integrable? Complete inte¬ 
grability requires n commitfmg Hamiltonians Hk : M —>■ R (i.e., their induced vector 
fields Xk = Xh^. commute, or equivalently, such that {Hi, Hj{ = 0 for the Poisson 
bracket induced by the 2-form w, {Hi, Hj{ = uj{XHi, Xh^) ) among themselves and 
with H itself. Finally, the Hamiltonians Hk should be functionally independent 
(i.e., their gradients should be linearly independent on a dense set of M). 

Very few dynamical systems are completely integrable, but these are especially 
important for being situations in which explicit computations may be performed. 
Indeed, one can make precise the idea that a generic Hamiltonian vector field does 
not have a second conserved quantity (recall that H itself is conserved along orbits). 
But this is just the opposite of what we expect from certain iterations in numerical 
analysis: to compute eigenvalues, for example, we expect to change something (an 
original matrix, an approximation of an eigenvector) without varying the objects 
being computed (the eigenvalues themselves). 

The Liouville-Arnold-Jost theorem states that, under appropriate hypothesis, 
the phase space of completely integrable equations foliates into invariant tori (i.e., 
products of lines and circles), given by levels of the conserved quantities Hk (fre¬ 
quently called the action variables in this context). In each torus, in angle variables, 
the evolution is just straight line motion. The angles vary smoothly at neighboring 
tori, and the global dynamics is mostly dependent on the arithmetic properties of 
the angular velocities. 

3.2 Toda flows in coadjoint orbits 

Using Flaschka’s change of variables, one might push forward the standard 2-form 
in R^" to iJo, the cone of Jacobi matrices with zero trace, thus converting it into 
a symplectic manifold. The surprising fact is that the resulting 2-form comes up 
from another construction of great interest, which we now describe briefly. 

A large class of symplectic manifolds is given by coadjoint orbits, equipped with 
the Lie-Kirillov-Poisson 2-form. Let G be a (finite dimensional) Lie group, 6 its Lie 
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algebra, and identify ©*, the dual of the Lie algebra, by means of a nondegenerate 
coupling 

(a, A) xe-)■ {a, A) € R. 

The coupling is bilinear and nondegeneracy means that the restrictions {a,.) and 
(., A), for a,A^Q, give rise to nonzero functionals respectively on ©* and ©. Thus, 
all linear functionals in © are of the form A {ct, A) for some a € ©*. 

For the Toda flow, start with G = U'^, the group of n x n upper triangular real 
matrices with positive diagonal entries. Then © = il is the vector space of real 
upper triangular matrices and we may identify ©* with S, the vector space of real, 
symmetric matrices, through the pairing (S,A) = trS'^l. 

The group G acts on itself by conjugation and on its Lie algebra Q by its deriva¬ 
tive at the origin, the adjoint action. For G = 17“'", the adjoint action is given by 
{g,U) G G X & Adg{U) = gUg~^ G ©. The coadjoint action from G x ©* to ©* 
satishes the natural duality, 

{Ad*g{a),A) = {a,Adg-i{A)), Va G ©*,Vy4 G ©. 

Back to the Toda context, 

Va G <S,VA G il(,7ld*a. A) = {a,Adg-i{A)) = iT:ag~^Ag = iT:{gag~^)A. 

The bad news is that gag~^ is not a symmetric matrix. Denote by sit the vector 
space of real, strictly upper triangular matrices. Clearly, for i? G sit and A G it we 
have trEA = 0. Consider the (unique) splitting 

M = S + E = UsAI + n^ixM, S' G 5, S G sit. 


We clearly have 

VA G it {Ad*a,A) = tT(gag~^)A = trn 5 (( 7 a( 7 “^)A, 

from which we finally obtain Ad*a = 115 ( 505 “^). 

By dehnition, M and have the same lower triangular part (and this in¬ 

cludes the diagonal). In particular, as observed by Adler [T] and Kostant [H], Jacobi 
matrices with fixed trace form a coadjoint orbit. 

We now recall the celebrated Lie-Kirillov-Poisson 2-form on a coadjoint orbit. 
Let Oa be the coadjoint orbit through a. Any vector a in the tangent space of Oa 
at a is the derivative at zero of Ad*^^f^^^^a fro some A in the Lie algebra ©. Let A 
and in © give rise to tangent vectors a and b. Set 

u}aia,b) = -{a, [A,B]). 

The minus sign is innocuous, but there is so much to prove here. First, it is not 
clear that w is well defined: other elements in © might give rise to the same tangent 
vectors at a. More, w has to be proven nondegenerate and closed. We will take all 
those issues for granted. 

Instead, we continue with the computations related to the Toda flow. In this 
case, the curve Ad*^^^^^s^a = II^ e*^ae~*^ has the tangent vector Eg [A, S'] at a. 
The 2-form is given by 

a;s(n 5 [A,S],n 5 [i 3 ,S]) = -trS[A,B]. 

A simple computation shows that the initial Hamiltonian H{x,y) for the Toda 
lattice in physical variables converts to H{S) = trS^/2 in Flaschka’s variables. 

Proposition 3 The Toda lattice is the vector field Xh associated to the Hamilto¬ 
nian H{S) = trS ^/2 defined on a coadjoint orbit Os- 
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Proof: We search for a vector field Xh = n5[17, S'] for which, given any vector 
field F = Ids [V, S] , we must have 

los{Xh, Y) =-tv S[U, F] = Y{H) = - tr SH^ [F, S]. 

Let 21 be the vector space (Lie algebra!) of real, skew symmetric matrices. We 
consider another splitting, 

M = A + U = U,kM + UupM, A G 2t, C/ e it 
and make use of the orthogonalities 21 _L S and it _L sit: 

tr S [U, F] = tr S Ug [F, S] = tr n„pS Ug [F, S] 

= trn„pS [F, S] = tr[ S, n„pS] F = - tr[ S, S], 

so that — tr[S, U] V = tr[ S, Ilsfc S] F, for all F G it, which is equivalent to say 
that [S, F] = [ S, Ilsfe S] + E, for some E G sit, so that 

Xh = n5[s, u] = Ug[ s, n,fc s] = [ s, n,fc s]. 

We are thus led to a very geometric explanation for the fact that the solution 
of the equation (J), J'{t) = [J{t),B{J{t))] = [J(f), Ilsfe is always a Jacobi 
matrix, if J(0) is: it simply does not leave the coadjoint orbit Oj(o)- 

A similar computation obtains the commutativity of the Hamiltonians Hk{J) = 
tr k = 1,... ,n). Indeed, the computation above generalizes to yield the Adler- 
Koslant-Symes theorem ([39]), a criterion to identify families of commuting Hamil¬ 
tonians. For the Toda lattice on Jacobi matrices with trace given by the initial 
condition J(0), the iJfe’s form a family of n — 1 commuting flows which are func¬ 
tional independent throughout Oj(o), since Jacobi matrices have simple spectrum 
from Lemma O Thus, the original Toda vector field and any other vector field 
given by a Hamiltonian of the form iJ/(J) = tr/(J) are completely integrable. 

From Theorem 12.21 the related invariant Liouville-Arnold tori, given by sets of 
n X n Jacobi matrices with fixed spectrum, is diffeomorphic to the set of possible 
choices for norming constants — the positive octant of the unit sphere in 

3.3 Toda flows in larger phase spaces 

From the computations above, the same Hamiltonian H = tr iS'^/2 induces equation 
{S) on larger coadjoint orbits and the proof of complete integrability on these phase 
spaces requires many more conserved quantities. The problem was considered in 
two cases. In |8], the authors consider the orbits of maximal dimension given by 
2[n^/4] in5. The new commuting Hamiltonians H{S) are obtained by chopping: 
they are the (symmetric functions of the) roots of the determinants of the matrices 
obtained by removing the first k rows and last k columns of the matrix S — XL In 
analogy to Moser’s computations in the Jacobi orbit, angle variables are essentially 
the first components of (generalized) eigenvectors. 

The generic phase space for real, nonsymmetric matrices is handled in [9]. First, 
one needs to interpret the Toda equation as a Hamiltonian on an appropriate coad¬ 
joint orbit of dimension — n of a Lie group given by a semidirect product. At a 
matrix M, the new commuting Hamiltonians are the coefficients of the polynomial 
p{z, A) = det(M — zM'^ — A), so a Riemann surface p{z, A) = 0 is associated to M 
and is invariant under the Toda flow. The existence of such additional structure 
has been known since the first studies of the periodic Toda flow, where the particles 
move on a circle instead of in the line ( [45] . |24] . |25)L For nonsymmetric matrices, 
the additional angle variables are obtained from an extension of the Abel-Jacobi 
map, by integrating a set of explicit meromorphic 1-forms on the surface along 
specific divisors, related again by generalized eigenvectors. 


3.4 The isospectral manifold of tridiagonal matrices 

The set Ji^ of Jacobi matrices with a given simple spectrum 

Ai < A 2 < ... < An, A = diag(Ai, A 2 ,..., A„) 

is diffeomorphic to There is a natural enlargement of this set — define 7a, 

the set of real, symmetric tridiagonal matrices with spectrum A. This set is actually 
a compact manifold ([43]). 





Figure 1: Ja, for n = 3 and 4 

The first step in the understanding of Ta is taking the closure (within the space 
of real,symmetric matrices) of Jt-^. It turns out that J7 a has an interesting combi¬ 
natorial structure, which we now describe. 

For matrix dimensions n = 3 and n = 4 (and any choices of different eigenvalues) 
this set is homeomorphic to the polytopes in the figure. The vertices correspond 
to the n! diagonal matrices with the same spectrum as A. For n = 3, the diagram 
consists of matrices with eigenvalues 2, 4 and 8. Vertices are diagonal matrices, in a 
self-explanatory notation. Edges correspond to matrices having a single zero in the 
main off-diagonal entries. Thus the top edge consists of matrices with (1,2) entry 
equal to zero; the eigenvalue 8 is trapped at entry (1,1) and the bottom 2x2 block 
consists of Jacobi matrices with eigenvalues 4 and 2, whose closure is homeomorphic 
to the whole edge. Edges are invariant under the Toda flow, and arrows indicate 
the sense of the flow. All interior points have the same a and uj limits. 

For n = 4 the boundary still consists of points with zero off-diagonal entries, 
which split the matrix in two kinds of sets, eight of which are homeomorphic to the 
3x3 case (hexagons), and six which are homeomorphic to the product of two 2x2 
blocks with fixed spectrum, the quadrilaterals. Again, the curves represent some 
orbits of the Toda flow. For the general case (Sa), i), define the permutohedron 
7^A CC)nV7r^5’n {(^77(1) 7 ■ ■ ■ ; ) } ■ 

Theorem 3.1 The map from J7a to the permutohedron Va 
J = Q^AQ ^ diagiQAQ^) € Va} 


is a homeomorphism. 

The existence of this homeomorphism in [33] led Bloch, Flaschka and Ratiu 
([6]) to search for a proof of the statement in terms of Atiyah’s theorem on the 
convexity of the image of Hamiltonian torus actions (|2]). Their proof starts with 
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(5,7.4) 


(7,5,4) 



Figure 2: Gluing four copies of Ja 


the identification of the homeomorphism as a momentum map, given by the formula 
above. A simple proof was later presented in [81) . 

From to the full isospectral manifold 7a, it’s a game of mirrors (or, more 
precisely, the construction of an appropriate Coxeter group ([13]))- In a nutshell, 
dropping the signs of the off-diagonal entries of a real, symmetric tridiagonal matrix 
(preserving its symmetry!) does not change the eigenvalues — this is something that 
numerical analysts use systematically: to compute eigenvalues of a matrix in 7a, 
it suffices to handle Jacobi matrices. In particular, the sets of matrices in 7a with 
nonzero entries split into 2"“^ connected components, all isomorphic to J7 a. To get 
7a, one takes the closure of these components and glues them along faces which 
are naturally identified. For n = 3 and A = {7,5,4}, Figure [2] shows some edges 
already glued. Edges along the boundary have to be identified: boundary vertices 
must be the same, together with a sign (which?). The resulting bitorus is drawn 
in Figure [3] so as to emphasize the boundaries of the hexagons (more about this 
picture on Section [6]). 



Figure 3: For n = 3, 7a is a bitorus 

The universal cover of 7a is and the homology groups are torsion free 

(031, HZ]). The sum of the partial traces, F{S) = — i + l)sii is a perfect 

Morse function: the equilibria are the diagonal matrices, and the Betti numbers are 
obtained by counting how many such matrices there are with a given signature. 
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4 A missed opportunity 

We present another trail ([3T]) leading to the Toda equation: it might have antici¬ 
pated the study of these equations by fifty years. 

We fix notation. Let A4(n) consist of the n x n real matrices, SO{n) C A4(n) 
be the set of orthogonal matrices with positive determinant and (n) C A4 (n) be 
the set of upper triangular matrices with positive diagonal entries. 

Clearly, both Lie groups act on M{n) by conjugation. For M G A4(n}, set 

02i = {Q'^MQ, Q e S'O(n)}, = {RMR-\ R G U+{n)} 

be respectively the orthogonal and upper triangular orbit through M. Clearly 
and are connected manifolds of matrices with the same eigenvalues as M. 

For convenience, let M have positive simple spectrum (i.e., positive distinct 
eigenvalues) — we are interested in Xm — Om 1^ 

The dimension count might suggest that Xm is trivial, but this is not the case. 
Indeed, for a matrix X G S, there are Q G SO{n) and R G U{n) such that 

X = Q'^MQ = RMR-^ => {QR)M = M{QR). 

In words, the matrix QR commutes with M. Since M has simple spectrum, 
QR is a function of M, in the sense that there exists / : ffi.—>■ K. — {0} such that 
QR = f{M). Notice that that Q and R are real invertible matrices. This is the 
QR factorization of f{M) (see Appendix 1 for details). In general, a real matrix 
X with positive determinant admits a unique factorization X = where 

[M]q G SO{n) and [M]n G U^{n) . 

Adding up, given / : R —>• R+, we obtain a matrix ^/)(/, M) G Xm, 

V>(/, M) = [ f{M) ]l M [ f{M) ]q = [ /(M) M [ f{M) . 

4.1 Commuting flows in Xm 

The notation ipi^f, M) was introduced to suggest that we are close to a flow on Xm- 

Proposition 4 Let f,g:^^ R.'*" and M be a real matrix with real, simple spec¬ 
trum. Then 

ipifg,M) = 'ip{gf,M) = il;{g,ip{f, M)). 

Proof: It suffices to show that 

[ f{M) g{M) ]q = [ fiM) ]q [ gm,M)) ]q. 

Now, from the uniqueness of the QR decomposition, for X with detX > 0, 

[ QX ]q = Q[X ]q and [ XR ]q = [ X ]q, for any Q G SO{n), R G iY+(n). 
Using that f{PXP~^) = Pf{X)P~^ and the first equality just above, 

[ gm,M)) ]q = [ g{[ fiM)) M [ fiM)) ]q) ]q = 

[ [ fiM) giM) [ fiM) ]q ] Q = [ fiM) [ giM) [ fiM) ]q ] 

Go back to the beginning and use the second equality: 

[ fiM) ]q [ gm,M)) ]q = [ giM) [ fiM) ]q = 

[ giM) [ fiM) ]q [ fiM) ]q = [ giM) /(M)]^ = [ (/g)(M)]Q. 

The underlying flow should now be evident: this is the content of the next result. 
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Corollary 1 Let M, /, Im = '"I V’ above. The map 

t, XgRxTm -.='il){e*^,X) 

is a globally defined flow on Im- Flows associated to functions f and g commute. 

4.2 Diferential equations for the flows 

Take derivatives of (j)f{t,M) to obtain the underlying differential equation. Set 
Q{t) = [ eFM) = [ e/(^) and M(0) = M. 

M{t) = (j>{t , M) = Q{tf M(0) Qit) = R{t) M{0) R-\t) {*) 

where = Q(^t) R{t), so that, dropping time dependence, 

M' = (Q^y M(0) Q + Q^M(O) Q' , 

/(M(0)) = /(M(0)) Q R =Q' R + Q R'. 

To get a differential equation, eliminate M(0) = Q{f)M{f){Q{t)Y'■. 

M'= {Q^)'QM+ MQ'^ Q’ and Q f{M{0)) R = Q'R +Q R'. 

Since Q{t) is a curve of orthogonal matrices, A{t) = Q' is a curve of skew 
symmetric matrices. Thus 

M' = -AM + MA = [M, A] and f{M) = Q' + R'R-^. 

We are almost done: the matrix f{M) in the last equation is a sum of a skew 
symmetric Q' and and upper triangular R'R~^, which bring to mind the factor¬ 
ization M = Ilsfe M -I- IlupM from Section [321 Here H^fc M and M have the same 
strictly lower triangular part, which determines Hsfe M. Adding up, 

M'{t) = [ M{t),A{t) ] = [ M(t), Usk f{M{t)) ] . 

The case f{x) = a; is the Toda lattice after Flaschka’s change of variables. Equation 
(*) is Symes’s solution by factorization to the differential equation ([39]). 

The fact that the Toda flow admits two different formulas by factorization, given 
by orthogonal and upper triangular conjugation of the initial condition, immediately 
implies that the evolution preserves the real tridiagonal symmetric form of the 
initial condition. Indeed, orthogonal conjugations preserve symmetry, and upper 
triangular conjugations preserve the upper Hessenberg form (i.e., the only nonzero 
entries below the diagonal lie in the subdiagonal of entries (i -I- l,f)). 

5 Toda, QR and other algorithms 

In the fifties, Francis m came up with the QR algorithm to compute eigenvalues 
of symmetric matrices. Say S' = S'o is a real, symmetric matrix of positive simple 
spectrum. Consider the alternation of QR decompositions and reorderings. 

So = QoRo, Si = RoQo = Qo SqQq = RqSqRq^. 

Clearly, Si € Osg, so So and Si are both symmetric with a common spectrum: 
the QR step is a diffeomorphism from Osq to itself. The remarkable thing about it 
is that iteration of this map converges to a diagonal matrix A — since a QR step 
clearly preserves spectrum, the diagonal entries of A are the eigenvalues of Sq! 
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It was Moser, again, who drew Deift’s attention to Symes’s beautiful connection 
between the Toda flow and the QR algorithm [ID]: at integer times n, the solution 
J{t) of the Toda differential equation satisfies exp J(n) = Em where En is the n-th 
term in the QR sequence starting from Eq = expJo. On the other hand, since the 
Toda equation is one within a family of flows of the kind J = [f,IlskewfiJ)]-, one 
may fudge with the functional parameter l|12jl and get a simpler relationship. 

Theorem 5.1 The flow associated to S = [/, Ilsfceu, In S'] interpolates the QR se¬ 
quence Sn, where So = S(0). 


Proof: From Section lT2l the solution of the equation is given by 

S{t) = Q*(t)S{0)Q{t), where exptlnS(O) = Q(t)R{t). 


Thus, So = J(0) = (5(1)S(1) and S(l) = (5*(1)S(0)(5(1) = Ji, the matrix obtained 
from So from a QR step. 

Since a Toda flow interpolates the QR iteration, the convergence properties of 
Toda flows are also satisfied by the iteration. The vocabulary of dynamical systems 
clarifies certain eigenvalue computations. The right side of Figure |T] represents Ja, 
the closure of the set of 4 x 4 Jacobi matrices with eigenvalues 1,2,3 and 4. The 
vertices correspond to diagonal matrices, which are equilibria for the Toda vector 
field. The vertex / is a source, A is a sink and the remaining vertices are saddles 
with different signatures (i.e., dimensions of the unstable manifold). The presence 
of saddles explains why orbits bifurcate close to some vertices (say, H, E and E) in 
the neighborhood of which an orbit spends a long time (i.e., many QR iterations), 
a fact that was known in the numerical literature as root disorder. 

Numerical analysts might have realized a long time ago that the QR algorithm 
is the integer evaluation of a flow. Indeed, it has been known for decades that one 
can obtain directly the matrix of the QR iteration starting with a symmetric 
matrix S'o is given by Sn = QnSoQn = RnSoRf^ where Qn and Rn are obtained 
from the QR factorization Sq = QnRn- Morally (and indeed correctly), the step 


Sl/n = Ql/nSoQl/n = Rl/nSoRf/n’ = Ql/nRl/n 


^1/n 


is an n-th root of the usual QR step (in the sense that n such steps yield the usual 
QR step). Now, to obtain the interpolating flow (which belongs to the Toda family, 
as we saw) simply compute 


lim n{Si/n - So) = lim n{Ql,.^SoQi/n - Sq)- (*) 

n—¥oo ' n^oo ' ' 

n n n 

. We are now ready to take the limit (*): the interpo- 
ln(a:), is 

S'it) = ln5]. 


Up to order I/n, 


In So 


and thus Qi/„ ~ / J- Il^fc 
lating flow, for which f{x) = 


5.1 Choleski and SVD 

There is nothing special about QR factorizations: other factorizations give rise to 
flows which are very similar to the Toda lattice. 

A Cholesky factorization M = MlMjj decomposes a matrix M as a product of 
lower and upper triangular matrices Ml and Mu with the same positive diagonal. 
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The factorization can be performed (uniquely) for matrices in having upper 
principal minors with positive determinant. Now (0), on G+ define the product 
g * h = hLghjj. For an appropriate coupling, the induced dual Lie algebra is the 
phase space which accommodates the Cholesky iteration, 

Alji — l/jiUrii ^n+1 — UjiL-fii ^ — 0, 1, . . . 

and its continuous interpolation (notice that blowups may happen). The Lie bracket 
associated to this group structure is an example of the so called R-matrix formalism 
applied to the standard matrix Lie bracket, but this is another story. 

Given a real matrix M, its singular values are the (nonnegative) lengths of the 
semi-axis of the ellipsoid obtained by applying M to the unit (Euclidean) sphere. 
A singular value decomposition of M is a product M = QYJJ, where Q and U 
are orthogonal matrices and E is a diagonal matrix having the singular values as 
diagonal entries. 

There is an efficient algorithm, similar to QR, to compute singular values of 
tridiagonal matrices, which was shown by Demmel and Kahan to have remarkable 
stability properties with respect to relative errors ([E])- In [7], these properties 
were studied from a symplectic setup: the appropriate phase space is chosen taking 
into account the specific concern with relative errors. The Jacobian M{i,j) of the 
map sending a matrix at step i to its (discrete) evolution at time j is analyzed 
using Krein’s perturbation theory for symplectic matrices. A number of properties 
arise, which are responsible for the good performance of the algorithm: for large 
i, the spectrum of M{i,j) is simple, lies in the unit circle and M(i,j) converges 
to a limit, explicitly computed. The agreement between theoretical estimates and 
experiments is remarkable: the rate between computed and estimated error was 
never larger than 8, independent of dimension. 

6 Bidiagonal coordinates 

Norming constants in Section 12.21 have a drawback: they do not cover the limit 
points of algorithms which converge to reduced matrices (i.e., tridiagonal matrices 
with some main off-diagonal entries equal to zero), like the Toda flows and QR type 
algorithms. This problem has been circumvented by the introduction of bidiagonal 
coordinates ([28]). As an extra bonus, bidiagonal coordinates provide an atlas for 
7a- The construction goes as follows. 

Let (n) denote the group of lower triangular matrices with unit diagonal en¬ 
tries. For M S Ai{n), the LU positive factorization, when it exists, is M = LU, 
where L G dl}{n) and U G Clearly, this happens if and only if the determi¬ 

nants of the upper principal minors of M are positive (more on the appendix). 

Let Sn be the symmetric group on the set {1, 2, ..., n}. For tt G S'„, define 

Att diag( ■ ■ ■ ; '^7r(n))- 

There is one chart of 7a for each permutation tt. Each chart has for domain 
the set C 7a consisting of matrices T = Q* A-n-Qjr for which there exists an 
orthogonal matrix admitting an LU positive factorization This is 

not as restrictive as it looks: if T = Q^At^Qtt, another spectral decomposition is 
given by r = Q’^E At^EQ^^, where A is a diagonal matrix with entries equal to ±1 
along the diagonal: one may use E to force the positivity of the determinants of 
the principal minors of EQ^^, provided that the corresponding determinants of Q-n- 
are nonzero. Notice that the request that Qt^ admits an LU positive factorization 
gives rise to a unique spectral decomposition T = Q* A^n-Q^r- We then have 

T = Ql A^Q^ = iLM~" A^L^U^ = U-^B^U^ 
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where 

= L-^ = U^TU-K 

Theorem 6.1 The matrix Bt^ is lower hidiagonal with diagonal The principal 
off-diagonal entries Pf = (i?-,r)fc+i,fc define a diffeomorphism : U\ 

Each domain hCff is an open, dense set IVf, containing one diagonal matrix. The 
charts : ITf —>■ K"“^,7r G Sn} form an atlas for Ta- The signs of Pf and ^ 
are equal and their quotient goes to one, when one of them goes to zero. 

Proof: From the expressions for it is simultaneously lower triangular and 

upper Hessenberg. So it is actually lower bidiagonal, with the same spectrum as 

Att- 

To show that the chart is a diffeomorphism to consider the construction of 

its inverse. Build Bt^ out of off-diagonal entries and (distinct) eigenvalues A^(i). 
Diagonalize Bt^ = and get out of the QR factorization = Q^^Rm 

so that, automatically, Qjr admits an LU positive factorization, hence G Uff. 
Finally set T = A^Q^r- 

Since U G W+(n), the equation = UT^TUf^ gives that the signs of Pf and 
Tk+i,k are equal. The remaining statements are left to the reader. 

Figure [2] is an example for n = 3: here A,r = diag(7,4, 5) and WJ is the interior 
of the polygon with boundary given by the unglued edges (the glued edges belong 
toUff). Figure [3] was obtained using such charts: the standard norming constants 
would distort too much the picture (and degenerate completely) at boundaries of 
(signed) Jacobi matrices. Once eigenvalues are fixed, the bitorus 7a lies in the 
intersection of a hyperplane (of matrices with the same trace as S) and a sphere 
(same sum of squared eigenvalues), and the figure is the image of a conformal 
projection of 7a in . 

Another remarkable property of bidiagonal coordinates is that their evolution 
under the Toda equations manages to be even simpler that the evolution of the 
standard inverse variables. We consider the Toda vector fields on real, symmetric, 
tridiagonal matrices, where the time dependence is explicit, 

r'(t) = [T(t),n,fe(/(T(t)))]. IT) 

Proposition 5 Fix tt, take T(0) G The chart domain IA\ is invariant under 
equation (T). The evolution of the bidiagonal coordinates is 

B^t) = [B.(t),-/(A^)], *.e., (Pint) = (/(A^(,+i)) -/(A^(q)) p^t) . 

Proof: We first prove that (T) leaves Ltfl invariant. Take T(0) G Ufl: omit¬ 

ting the permutation tt, we have T(0) = QqAQq, where Qo has an LU positive 
decomposition. Solve (T) as in Sectionfor exp(t/(T(0))) = Q(t)R{t), 

T{t) = Q{f)* T(0) Q(t) = Q{t)* Q*o A Qo Q{t) 

We have to prove that Q{t) = Qq Q{t) admits an LU positive factorization: 

Q(t) = Qo exp(t/(T(0))) R~^{t) = Qo exp(t/(T(0))) Qq Qo R~^it) 

= Qo exp(t/(Qo AQo))Qo QoR~^{t) = exp(t/(A)) Qo i?"^(t). 

The upper principal minors of Q{t) and Qo have the same signs, since exp(t/(A)) 
is a positive diagonal matrix and R~^{t) G ld^{n). 

We now consider bidiagonal coordinates. Clearly 

T{t) = Q*{t)Q*AQoQ{t) = Q*{t) AQ{t), Q* If) Q'it) = Q* {t)Q'(t). 
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Now, B{t) = L~^(t) AL{t) and we obtain B'{t) = [B{t),L~^(t) L'{t)\ by now 
familiar computations. The matrix L is obtained by the LU positive factorization 
Q{t) = L{t) U{t), so that, dropping the time dependence, Q' = L'U + LU’ yields 

L-^Q’U-'^ = L-^L' + U'U-^. 

Consider the split M = IIsj M + n„M of a matrix M into strictly lower and upper 
triangular parts. Since L has a diagonal of ones, L~^L' is strictly lower triangular. 
The solution by factorization in Section gives Q* Q' = Ilsfc so that 

L-^L' = L-^Q’u-^ = u{Usk fiT)) u-\ 

Now, the matrices M and 11^^ M have the same strictly lower triangular part and 
U is upper triangular, so 

L-^L' = uf{T)u-^ = f{UTU-^) = n,z /(s). 

Adding up, 

B' = [B, L-^L’] = [B, f{B)] = [B, f{B) - /(A)] = [B, -/(A)] 


The proof requires interpretation in the case f{x) = In a;, which is of relevance 
for QR interpolation: we need exp(t/(A)) to be a positive diagonal matrix — in 
this case, the diagonal entries must be equal to |Afe|. 

Bidiagonal coordinates are especially convenient to study asymptotic behavior 
of Toda flows. As an application, the reader may find in |28] a rather natural 
computation of their scattering map, described in Section 12.21 by filling up the 
following inevitable outline. Recall from Figure 1 that Toda flows starting from 
Jacobi matrices have for w and a limits the diagonal matrices associated to (i) = i 
and TTaii) = n — i + 1. On each chart, the Toda evolution in bidiagonal coordinates 
is simple. The change of charts required to keep track of both extremes of an orbit 
is equally simple. 

An appropriate extension of this formalism provides charts on isospectral man¬ 
ifolds of real and complex matrices with given profile^ a natural extension of the 
concept of tridiagonality — a text is under preparation (mi)- 


6.1 QR steps with shifts 

Numerical analysts have ways to speed up the convergence of QR. The original QR 
step is just the choice f{x) = Inx in the family 

Sn i-t Sn+i = [exp/(S'o)]Q S'o [exp/(S'o)]Q = [exp/(S'o)]fl S'o [exp/(S'o)]^^ 

Bidiagonal coordinates may be used to indicate interesting alternatives: interpret a 
step as a time 1 map for a differential equation and integrate the trivial flow which 
describe evolutions in Proposition [S] For f{x) = hig(x), taking into account the 
caveat after its proof, the change of the bidiagonal coordinates of a matrix T G 
under a step is 


(/9 


TT 
1 5 


>Pn-l) 


g(^7r(2)) 

5('^7r(l)) 


PI, 


9{^7r{n )) 
g(A7r(n—1)) 



We are interested in functions g for which some becomes small — in this case, 
the corresponding matrix essentially decouples in two smaller tridiagonal matrices, 
for which the computation of eigenvalues is simpler. Taking into account the de¬ 
nominators in the formula, a natural possibility is a function g which equals zero 
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at A^(„) but does not vanish at the other eigenvalues. Keep in mind that, unfortu¬ 
nately, we do not know the eigenvalues. In particular, it is especially hard to obtain 
such a g which would reduce drastically a centrally located (i.e., k nj2). 

Notice that, from Theorem 16.11 the matrix T is Jacobi if and only if its bidiag¬ 
onal coordinates are positive, and the sign is preserved under general QR steps. 
On the other hand, removing the absolute values in the formula shows that one can 
replace iterations lying within Jacobi matrices to iterations on tridiagonal matrices 
without changing their asymptotic properties: only the off-diagonal entries eventu¬ 
ally have different signs. The new iteration is smooth on the isospectral manifold 
7a, and convergence rates may be obtained using Taylor expansions. 

Typically, along an iteration of an algorithm searching for eigenvalues, one has 
good approximations s for one of them: a natural choice is g(x) = x — s, the iteration 
with shift s. The bidiagonal coordinates, in this case, change as follows: 



Numerical analysts frequently do not wait for a good approximation s. There 
are different shift strategies ( [34] 1 — we consider the Rayleigh quotient shift for 
which s = Tn^ni and the Wilkinson shift: compute the eigenvalues of T, the bottom 
2x2 principal minor of T, and take for s the one which is closer to Tn^n- 

Under the Rayleigh shift, once iteration approaches convergence, a simple Tay¬ 
lor expansion shows that the bottom entry b = /3^_i converges cubically to zero, in 
the sense that 6™+^ = 0((6™)^) for bottom entries at consecutive iterations. Once 
b is small enough, the matrix undergoes deflation: is declared a good approxi¬ 

mation of an eigenvalue and the algorithm proceeds with the top (n — 1) x (n — 1) 
block. There is one catch however: for some initial conditions, the Rayleigh shift 
gives rise to periodic orbits. 

This does not happen for the Wilkinson shift. The dynamics in this case is 
richer: we describe the results but may only indicate [29] and [30] for proofs. 

Theorem 6.2 For a generic initial condition, a Wilkinson iteration leads to cubic 
convergence ofb = If the original matrix T has no three eigenvalues 

in arithmetic progression, this is always the case. Otherwise, there may he a Cantor¬ 
like set of initial conditions for which iteration is quadratic. 

In a nutshell, this peculiar behavior is caused by the discontinuities in the def¬ 
inition of s, at points where is equidistant from the eigenvalues of the block 
T. The Cantor-like set in the statement of the theorem consists of matrices all of 
whose iterations lie in this situation. Typically, this never happens or happens just 
at few steps and convergence is cubic. 

7 Lax pairs beyond Toda 

Forty years of contributions from a large community greatly increased our under¬ 
standing of the Toda flow, which sometimes is undistinguishable from the more 
general Lax pair evolution. Mutations are abundant: for a beautiful starting point 
about different representations of the Toda flow (and other geometric issues which 
may studied through them), the reader should refer to [23] . 

In this text, we emphasized examples over theory. In order to indicate the 
versatility of the concepts which have been presented, we close with a final example. 
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7.1 The billiard on an ellipsoid 

Moser and Veselov [33] introduced a very interesting formulation of the billiard 
on an ellipsoid E = {x : {x,C~^x) < 1}. We fix notation: a ball moves along a 
direction j/o (a unit vector), hits the boundary dE at xq, and leaves in direction yi 
until it hits Xi. They presented the billiard map 'I'(a;o, j/o) = (a^i, yi) in the following 
fashion. Write the (matrix) polynomial 

To(A) = yo ® yo + Axo A yo - A^C^, 


factor in linear terms, 

To(A) = (AC + yo (8) ^o){-XC + ® yo), Co = C~^xo, ||Co|| = 1, 

intertwine factors as in the QR algorithm and factor again, 

'^o(A) = (—AC + Co <8 yo){XC + yo ® Co) = (AC + yo <8 Co)(~AC + Co ® yo)) 
with y'g = V'o, IIColl = 1- Set C<(a;o,yo) = (-Cyi, C'^Xo): then = -(fp. 

In dni , the analogy of the Moser-Veselov formalism to the QR iteration — better 
still, to the Cholesky iteration — is taken literally. The relevant group now is 
C+, consisting of loops 7 : iR —>• GL{n,C) which are smooth at 00, contract to 
the identity loop, are positive diagonal at 00, satisfy the reality condition 7(A) = 
7 (A) and admit a (unique) Riemann-Hilbert factorization 7 = 7 l 7 k, where ■jl 
and 7 i{ have analytic continuations to the left and right side of the imaginary axis 
respectively, and same diagonal values at 00. The group operation is 7 * <5 = SljSr. 
Take an initial data for the billiard evolution and construct the rational matrix 
function Ao(A) = Lo(A)/(l — A^). 

One has to circumvent a technical difficulty: Ao(0) is not invertible — it is a rank 
one matrix! Still, as in the finite dimensional cases, the evolution is interpolated by 
a differential equation. This time, the solution formula by factorization involves a 
Riemann-Hilbert problem, with a mild singularity at A = 0. 


Appendix 1: The QR and LU factorizations 

Let M an invertible, real matrix. Then there is a unique Q G SO{n) and R € Ll{n) 
for which M = QR, the QR factorization of M. Indeed, let e;, i = 1,..., n be the 
canonical vectors: the reader should have no difficulty in showing that the subspaces 
generated by Mci and Qci, for i = 1 ,..., fc {k arbitrary) should be the same, if 
such a decomposition exists. Since the columns of Q are orthonormal, they must 
be vectors obtained by applying the Gram-Schmidt orthogonalization procedure to 
the columns of M sequentially. The entries of R are the coefficients used in the 
representation of Mci in terms of the columns Qei, ..., Qci. Since M is invertible, 
the process is feasible, and appropriate normalizations give rise to positive diagonal 
entries of R. 

Similarly, the (unique) LU decomposition of M is M = LU, where L £ E^{n) 
and U £ Li^{n). This decomposition exists if and only if the principal diagonal 
minors of M (i.e., the determinants of the k xk submatrices with entries in the 
intersection of the first k rows and columns, k = 1,... ,n) are strictly positive. In 
a nutshell, from Mk = L^Uk, it is clear that an inductive construction is at hand: 
the details are left to the reader or in standard texts ([2], [44] 1. 
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